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ABSTRACT 

We analyse a hydrodynamical simulation model for the recurrent heating of the central 
intracluster medium (ICM) by active galactic nuclei (AGN). Besides the self-gravity 
of the dark matter and gas components, our approach includes the radiative cooling 
and photoheating of the gas, as well as a subresolution multiphase model for star 
formation and supernova feedback. Additionally, we incorporate a periodic heating 
mechanism in the form of hot, buoyant bubbles, injected into the intragalactic medium 
(IGM) during the active phases of the accreting central AGN. We use simulations 
of isolated cluster halos of different masses to study the bubble dynamics and the 
heat transport into the IGM. We also apply our model to self-consistent cosmological 
simulations of the formation of galaxy clusters with a range of masses. Our numerical 
schemes explore a variety of different assumptions for the spatial configuration of 
AGN-driven bubbles, for their duty cycles and for the energy injection mechanism, in 
order to obtain better constraints on the underlying physical picture. We argue that 
AGN heating can substantially affect the properties of both the stellar and gaseous 
components of clusters of galaxies. Most importantly, it alters the properties of the 
central dominant (cD) galaxy by reducing the mass deposition rate of freshly cooled 
gas out of the ICM, thereby offering an energetically plausible solution to the cooling 
flow problem. At the same time, this leads to reduced or eliminated star formation in 
the central cD galaxy, giving it red stellar colours as observed. 

Key words: methods: numerical - galaxies: clusters: general - cooling flows - cos- 
mology: theory. 



1 INTRODUCTION 

Clusters of galaxies are the largest virialised objects in the 
Universe and are thought to contain a representative fraction 
of baryons (White et al., 1993). Most of these baryons can be 
found in the diffuse gas of the intracluster medium (ICM), 
which is directly observable in X-rays, making clusters of 
galaxies an almost ideal laboratory for studying the physi- 
cal processes that shape galaxies and halos in the Universe. 
Clusters of galaxies are also a useful cosmological probe (for 
a recent review see Voit, 2004), and therefore have been a 
prime target for theoretical modelling early on, both numer- 
ically and analytically. 

A first order approximation for the ICM is to represent 
it as an ideal, non-radiative gas. This leads to the predica- 
tions of scale invariant relations between X-ray luminosity, 
mass and temperature (Kaiser, 1986). However, it has long 
been established that the observed relations do not agree 
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in detail with these assumptions, e.g. the observed Ly^-T 
relation is much steeper than expected based on this sim- 
ple model. In addition, recent observations with radio and 
X-ray telescopes have revealed a stunning complexity of the 
ICM physics, including phenomena such as cold fronts, radio 
ghosts, cluster turbulence, and apparently nearly uniform 
enrichment to high metallicity. 

However, possibly the most puzzling observational fact 
is the "cooling-flow" problem. Since the cooling time in the 
central regions of galaxy clusters is smaller than the age 
of the clusters themselves, a central inflow of cool gas is 
expected to occur (e.g. Fabian & Nulsen, 1977; Cowie & 
Binney, 1977; Fabian, 1994). The rate of gas cooling can be 
estimated by energetic arguments if one assumes that X-ray 
cooling radiation is fed by the thermal reservoir of the hot 
cluster plasma. Based on this, the estimated rate of accre- 
tion onto the central galaxy is rather high in many cases (e.g. 
Fabian et al, 1984; White et al., 1994, 1997; Allen, 2000; 
Allen et al., 2001a), implying that a significant amount of gas 
cooler than 1 — 2 keV should be present in the centre. How- 
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ever, up to the present time, optical and X-ray observations 
have failed to detect the required amount of this cool gas, 
suggesting that it is simply not there (e.g. McNamara et al., 
2000; Peterson et al., 2001; Tamura et al., 2001; Balogh 
et al, 2001; Kaastra et al., 2001; Edge, 2001; Edge et al., 
2002; Edge & Frayer, 2003; Fabian et al., 2001; Bohringer 
et al, 2002; Peterson et al., 2003; Salome & Combes, 2004). 
The low current star formation rates of central galaxies (e.g. 
O'Connell & McNamara, 1989; Johnstone et al., 1987; Allen 
et al., 1995) provide additional support for the absence of 
strong cooling flows. Apparently, there must be a physical 
process that offsets the radiative cooling in the centre, pre- 
venting the gas from falling out of the ICM in a cooling 
flow. 

Theoretical studies have therefore often invoked some 
sort of non-gravitational heating to explain the cluster scal- 
ing relations (e.g. Kaiser, 1991; Navarro et al., 1995; Bower, 
1997; Tozzi & Norman, 2001; Borgani et al., 2001; Voit & 
Bryan, 2001; Babul et al., 2002; Voit et al., 2002, 2003; 
Oh & Benson, 2003; Tornatore et al., 2003; Borgani et al., 
2004). The main unsolved issue in these models remains 
the origin and nature of the physical sources that cause 
the extra-heating of the ICM. Perhaps the most obvious 
heat source is supernovae associated with star formation, 
but it seems questionable that they are able to supply the re- 
quired amount of feedback energy. Curiously, radiative cool- 
ing alone may also account for the steepness of the Ly^-T 
relation by eliminating gas more efficiently in low-mass sys- 
tems (e.g. Lewis et al., 2000; Voit & Bryan, 2001; Muan- 
wong et al, 2001; Yoshida et al., 2002; Wu & Xue, 2002; 
Voit et al., 2002; McCarthy et al., 2004), but this produces 
a drastic overprediction of the amount of cold gas (apart 
from a problem with the Lx~T zero-point) and is therefore 
disfavoured. Models that self-consistently incorporate SNe 
heating and radiative cooling processes are also only found 
to have limited success (e.g. Borgani et al., 2004). The over- 
cooling problem therefore remains unsolved. Another prob- 
lem is posed by simulated temperature profiles, which typi- 
cally exhibit a trend to increase towards the cluster centre, in 
disagreement with observational inferences (e.g. Allen et al., 
2001b; De Grandi & Molendi, 2002). 

An array of different physical hypothesis have been pro- 
posed to solve the cooling-flow paradox, including thermal 
conduction, magnetic fields, cosmic rays, and hot buoyant 
bubbles from AGN jets. Thermal conduction may in prin- 
ciple offset central cooling losses by inducing a heat cur- 
rent from outer, hotter regions of clusters (e.g. Narayan & 
Medvedev, 2001), provided the conductivity is not strongly 
suppressed by tangled magnetic fields. Analysis of static 
cluster models with conduction have been able to provide 
good matches to observed temperature profiles in some cases 
(e.g Voigt et al., 2002; Zakamska & Narayan, 2003; Voigt & 
Fabian, 2003) but detailed self-consistent numerical simula- 
tions which followed conduction still encountered the cooling 
flow problem (e.g. Jubelgas et al., 2004; Dolag et al., 2004), 
making it questionable whether this can be the real solution. 

The more widely favoured hypothesis is instead that the 
central AGN may supply the required amount of energy. Ac- 
cretion onto supermassive black holes is thought to liberate 
of order ~ 10% of the accreted rest mass energy, implying 
that even for low accretion rates onto a supermassive black 
hole, offsetting the cooling flows is energetically quite possi- 



ble. In fact, such accretion powers high-redshift quasars, the 
most luminous sources in the universe. Quasar activity is 
likely to be triggered by mergers of galaxies, where cold gas 
is forced to the nuclei by gravitational tidal forces. This ac- 
cretion and the associated quasar feedback has recently been 
incorporated into simulations, and shown to play a poten- 
tially important role in shaping the properties of elliptical 
galaxies (Springel et al., 2005a). 

In clusters of galaxies, however, it seems clear that the 
central AGN activity that causes bubbles is of a different 
nature, and needs not be triggered by galaxy mergers. Ob- 
servationally, many clusters of galaxies show evidence for 
X-ray cavities filled with radio plasma (e.g. Owen et al., 
2000; Blanton et al., 2001), which are thought to be inflated 
by relativist ic jets from the AGN. Theoretically, it has been 
shown that these bubbles may rise buoyantly and raise some 
of the central cool gas (e.g. Churazov et al., 2001), allowing 
it to mix with the hotter gas further out. Together with the 
accompanying mechanical and possibly viscous heating, this 
can then constitute an efficient feedback mechanism. 

In this paper, we focus on the phenomenology of this 
bubble feedback, without addressing the small scale physics 
of the accretion onto the black hole. This extends earlier 
simulation studies which all employed hydrodynamical mesh 
codes, but which focused exclusively on highly idealised 
cluster models (e.g. Churazov et al., 2001; Quilis et al., 
2001; Ruszkowski & Begelman, 2002; Churazov et al., 2002; 
Briiggen et al., 2002; Briiggen & Kaiser, 2002; Briiggen, 
2003; Nulsen et al., 2003; Dalla Vecchia et al., 2004; Hoeft 
& Briiggen, 2004). A first goal of our work is to demonstrate 
that such simulations are also possible with the smoothed 
particle hydrodynamics (SPH) technique, and give results 
consistent with earlier studies. This is important because 
the Lagrangian nature of SPH is ideal for cosmological sim- 
ulations of structure formation, and if applicable for bubble 
feedback, will allow us to carry out the first self-consistent 
cosmological simulations with AGN-driven bubble heating. 
An equally important goal of our simulations is to gain new 
insights into the efficiency of bubble feedback associated 
with AGN for modifying the thermodynamic state of the 
ICM and the properties of cluster galaxies over the course 
of cosmic history. Our modelling can hence inform semi- 
analytic models of galaxy formation that have just begun to 
include AGN feedback (e.g. Croton et al., 2005), and pro- 
vide crucial input for future hydrodynamic simulations that 
try to incorporate the growth of supermassive black holes 
both from the quasar- and the radio-mode. 

The outline of this paper is as follows. In Section 2, we 
describe the characteristics of our simulation code and the 
numerical method adopted to introduce bubble heating. In 
Section 3 we analyse the AGN heating in isolated galaxy 
halos, spanning a wide range in mass, and we present some 
Chandra- like photon images of simulated bubbles. The ef- 
fects of AGN heating in cosmological simulations of galaxy 
cluster formation is discussed in Section 4. Finally, in Sec- 
tion 5 we discuss successes and limitations of our model, and 
we present our conclusions. 
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2 METHODOLOGY 
2.1 Basic code properties 

Our simulations have been performed with the parallel 
TreeSPH-code GADGET- 2 (Springel, 2005; Springel et ah, 
2001b). We use the 'entropy formulation' for SPH suggested 
by Springel & Hernquist (2002), which manifestly conserves 
both energy and entropy when adaptive smoothing lengths 
are used. Besides gravitational and hydrodynamical pro- 
cesses, we include radiative cooling of the gas component, 
together with heating by a spatially uniform, time depen- 
dent UV background modelled as in Katz et al. (1996). The 
gas consists of an optically thin primordial plasma of hy- 
drogen and helium. In addition, a multiphase subresolution 
model for the treatment of star formation and associated 
feedback mechanisms has been adopted (Springel & Hern- 
quist, 2003). In this model, stars form from dense cold gas 
clouds assumed to reside at pressure equilibrium in a sur- 
rounding hot phase of the interstellar medium. Supernova 
explosions heat the hot medium and evaporate cold clouds, 
thereby providing a self-regulation cycle for star formation, 
and a net pressurisation for the highly overdense ISM. Addi- 
tionally, we use a simple prescription for metal enrichment, 
which assumes that each star-forming gas element can be 
locally approximated by a closed box model in which per- 
fect and instantaneous mixing of metals between cold clouds 
and ambient gas occurs, as explained in detail in Springel & 
Hernquist (2003). 



2.2 Phenomenological description of AGN 
heating in clusters 

Besides considering the physical processes already imple- 
mented in GADGET-2, we have implemented for this study 
a new model that accounts for heating by the central AGN 
in clusters of galaxies. This model does not attempt to pro- 
vide a fully self-consistent ab initio treatment of the com- 
plex physical processes related to accretion onto supermas- 
sive black holes in clusters and the associated AGN activity. 
Rather, we try to mimic the observed phenomenology of hot 
bubbles in clusters directly in our simulations, without ad- 
dressing the jet physics that presumably inflates the bubbles 
in the first place. We therefore assume as a starting point 
that such bubbles are generated during phases in which an 
AGN is "switched on" , and introduce them into the IGM in 
a phenomenological fashion. This allows us to study how the 
bubbles affect the properties of the central ICM as a func- 
tion of their characteristics, in particular with respect to 
distributing their energy content to the surrounding cooler 
gas. 

For definiteness, we assume in our model that a certain 
amount of thermal energy is injected in the form of centrally 
concentrated bubbles spaced in uniform time intervals. We 
parameterise this scheme in terms of the AGN duty cycle, 
the amount of energy Ehuh injected, and by the radius Rbub 
and distance dbub of the buoyant bubbles from the cluster 
centre, respectively. 

We first test our scheme for AGN-heating on isolated, 
axisymmetric halo models. These systems are clean labo- 
ratories which permit us to compare directly with analo- 
gous modelling in the literature (e.g. Churazov et al., 2001; 



Quilis et al., 2001; Dalla Vecchia et al., 2004), and hence 
to evaluate whether SPH is suitable for such simulations. 
Moreover, these simplified models give us the possibility to 
explore straightforwardly and with comparatively low com- 
putational cost a large number of cases. In this way we can 
investigate the importance of different physical parameters 
of the bubbles, thus constraining their dynamical evolution 
and the heat transport into the ICM. 

As a second step, we apply the model for bubble 
heating to fully self-consistent cosmological simulations of 
galaxy cluster formation. Here, we also investigate different 
redshift-dependent energy injection schemes, allowing us to 
gain some insight in how the AGN activity influences the 
hierarchical galaxy cluster growth and the characteristics of 
the central cluster galaxy, and to elucidate the relative im- 
portance of AGN heating with respect to the other physics 
included. We consider a set a galaxy clusters spanning a 
range in mass because we expect the efficiency of bubble 
heating to have a significant mass dependence. 

Both for isolated halos and in cosmological simulations, 
we explored two different schemes for spatially placing the 
bubbles around the cluster centres. In the first scheme, the 
bubbles are introduced randomly within a sphere with a 
radius given by dbub around the centre, while in the second 
approach, two symmetric bubbles are placed randomly along 
a fixed axis of length 2 x dbub, which has an orientation pre- 
served in time during subsequent bubble events. The latter 
hence mimics a situation where the AGN jet that inflates 
the bubbles has directional stability over time, which could 
arise due to some coupling with the host galaxy's angular 
momentum, for example. At the present time there is no 
clear evidence either way concerning what is the preferred 
scenario, therefore our main aim is to investigate the pos- 
sible differences in the ICM properties between these two 
bracketing scenarios. 

2.3 Constraining the model parameters 

Our choice for the values of Rbub and cfcub has been guided 
by observational constraints on X-ray cavities in clusters, 
and also by the values typically adopted in previous numeri- 
cal works, for easier comparison. For simplicity, we restricted 
most of our simulations to the case where the values of it^ub 
and dbub depend only on the mass of the halo under con- 
sideration, and on the redshift in the case of cosmological 
simulations. Specifically, we adopted 

l + z 



R huh oc M 20 o(z) 



1/3 



(Qoma + ^ + ^OA) 1 / 3 ' 



(1) 



where M2oo(z) is the virial mass of the host galaxy cluster 
at given redshift of AGN activity, and the same scaling has 
been adopted for dbub- For the simulations of isolated halos, 
we used the same dependence of Rbub and dbub on cluster 
mass, setting z = 0. 

We study multiple bubble injection events in order to 
analyze how AGN heating couples with radiative cooling 
losses over a sufficiently long time interval. Thus, our model- 
ing requires prescriptions both for the AGN duty cycle and 
for the time evolution of the energy content stored in the 
bubbles. However, most of the observed AGN-driven bub- 
bles are found at low redshifts (e.g. Birzan et al., 2004), and 
only recently some observational evidence for X-ray cavities 
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in more distant galaxy clusters has been found (McNamara 
et al., 2005). Therefore, the properties and presence of radio- 
bubbles at higher redshifts, and their evolution with time, 
are observationally rather unconstrained. We hence limit 
ourselves in this work to simple parametric prescriptions for 
the evolution of i^bub, derived from basic theoretical consid- 
erations and empirical laws, which hopefully bracket reality. 
Typically, we started injecting bubbles at redshift z = 3, 
which is the epoch that roughly corresponds to the peak of 
the comoving quasar space density, but we also tested an 
earlier epoch given by z = 6 for the start of the bubble ac- 
tivity. For our modeling of the evolution of i^bub with time, 
we adopted two scenarios with rather different behaviour. In 
the first one, most of the energy is released at late epochs, 
while in the second one, the bubble energy is coupled more 
closely to an assumed BH accretion rate (BHAR) model for 
the growth of the black hole population as a whole, such that 
the energy release is more pronounced at high redshifts. 

More specifically, our first model is loosely motivated 
by the Magorrian relationship, which implies Mbh oc <t 4 . 
A relation between the bubble mechanical luminosity and 
the black hole accretion rate, Mbh, can be derived by as- 
suming that only a small fraction of the total bolometric 
luminosity thermally couples with the ICM. Hence, Lbub = 
/ x Lboi = / x cMbhc 2 . The factor / sets the efficiency of 
thermal coupling with the ICM, and is typically assumed 
to lie in the range of 1-5%, while e is the radiative effi- 
ciency factor. Assuming that the mechanical luminosity for 
Eddington-limited accretion is directly proportional to the 
black hole mass, the energy content i^bub of the bubbles is 
then proportional to M200 (z) 4//3 , provided the mass of the 
central cluster galaxy scales self-similarly with the cluster 
mass. Hence, it follows that in this model the bubble en- 
ergy content is determined by the mass assembly of the host 
galaxy cluster with time. 

In our second scenario, we instead relate the amount of 
bubble energy to the average growth rate of supermassive 
central black holes. To describe the latter, we employ an 
estimate of the BHAR by Di Matteo et al. (2003), who give 
an analytic fit 



6exp[a(z - Zm)] 
b — a + aexp[6(z — z 



0] 



(2) 



for their numerical results, with the parameters a = 5/4, 
b = 3/2, z m = 4.8, and e BH = 3 x 10" 4 M o yr" 1 Mpc" 3 . 
Thus, for every duty cycle of AGN activity we can directly 
relate Mbh with i^bub in the following way, 



^bub - w w 2 

= / X € X C 



p(z) dz, 



(3) 



where a normalisation factor, E norm , has been introduced 
which we set such that the total energy injected over all 
duty cycles is the same in our two schemes. We note that 
the different temporal evolution of the BH mass in this ap- 
proach implies a significantly reduced energy content of the 
bubbles at low redshifts. Finally, there is still one free con- 
stant of integration which we choose by requiring that the 
assumed mass of the black hole is the same at z — in both 
scenarious. 

A number of observational and theoretical works (e.g. 
Birzan et al., 2004; Sanderson et al., 2005; McNamara et al., 
2005; Nulsen et al., 2005a,b; Donahue et al., 2005; Voit & 



Donahue, 2005) have constrained the plausible time interval 
between two successive bubble injection episodes to be of 
order of Atbub ~ 10 8 yrs. Clearly, Atbub could vary both for 
clusters of different mass and also in time, especially if the 
bubble activity is triggered by a self-regulated mechanism 
that operates between AGN feedback and the cooling flow. 
Nevertheless, given our simple phenomenological approach 
and lack of any better observational constraints, we adopt 
the value of Atbub = 10 8 yrs for all of our cluster simula- 
tions, independent of the cosmological epoch. 

While some of our prescriptions for bubble parameters 
are motivated by "quasar- like" phenomena, our models are 
really meant to reflect a mode of feedback by supermassive 
black holes different from that of ordinary quasars. Instead 
of being triggered by mergers and being fueled with dense 
and cold ISM gas, the bubbles are a model for the radio 
activity observed in clusters. Note that there are also newly 
emerging theoretical models (e.g. Croton et al., 2005; Chu- 
razov et al., 2005) on how both quasar activity at higher 
redshifts and AGN-driven radio bubbles at lower redshifts 
can be described within a common unified framework. We 
will discuss this possibility in more detail in our conclusions. 



3 AGN HEATING OF ISOLATED GALAXY 
CLUSTERS 

We here analyse simulations of isolated halos, consisting of 
a static NFW dark matter halo (Navarro et al., 1996, 1997) 
with a gaseous component that is initially in hydrostatic 
equilibrium and chosen to follow a density distribution sim- 
ilar to the NFW dark matter profile, but slightly softened 
at the centre according to 

fb ^0 Pcrit 



PgiT) 



(4) 



(r + r )/r s (l + r/r s ) 2 ' 

where ro is a parameter introduced to mimic a gas core 
radius. The baryonic fraction is given by fb, while p cr it is 
the critical density, £0 is the characteristic overdensity and 
r s is the scale radius. The gas follows the equation of state 
of an ideal monoatomic gas with adiabatic index 7 = 5/3. 
Besides, a certain amount of angular momentum has been 
imposed that can be quantified by the dimensionless spin 
parameter of a halo, 

A - J|£|V2 (5) 

where J represents the angular momentum, M is the halo 
mass, and E its total energy. 

The boundary conditions were chosen to be vacuum, 
i.e. both density and pressure are initially zero outside the 
virial radius (defined here as the radius enclosing a mean 
density equal to 200p C rit). We have simulated halos with a 
wide range of masses, with virial radii and concentration 
parameters as listed in Table 1. The baryonic fraction, fb = 
0.12, the spin parameter, A = 0.05, and ro = 0.02i?2oo were 
kept fixed for all the halos. When evolved without radiative 
cooling, these initial models are perfectly stable for more 
than 1/4 of the Hubble time, as we explicitly checked. This 
is the timespan we subsequently consider in all our non- 
trivial simulations, both for the case with cooling and star 
formation only, and also for the case with additional AGN 
heating. 
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M 20 o [h- 1 M Q ] 


#200 [h 1 kpc] 


c 


Ngas 


m g as [h 1 M Q ] 


e [fc-ikpc] 


10 12 


206 


12.0 


3 x 10 5 


4.0 x 10 5 


1.0 


10 13 


444 


6.5 


3 x 10 5 


4.0 x 10 6 


2.0 


10 14 


957 


8.0 


3 x 10 5 


4.0 x 10 7 


5.0 


10 15 


2063 


5.0 


3 x 10 5 


4.0 x 10 8 


10.0 


10 15 


2063 


5.0 


1 x 10 6 


1.2 x 10 8 


6.5 



Table 1. Numerical parameters of the isolated galaxy clusters. The first two columns give the virial mass and radius of the halos, 
evaluated at 200p cr it- The assumed values for the concentration parameter are in the third column, while the initial number and the 
mass of the gas particles is shown in the fourth and the fifth columns, respectively. The mass of the star particles is half that of the 
gas particles, because we set the number of generations of star particles that a gas particle may produce to two. Note that there are no 
parameters for the dark matter particles in these run, because we modelled the dark halo with a static NFW potential. Finally, in the 
last column, the gravitational softening length e for the gas and star particles is given. 



For the 1O 15 /i _1 M isolated cluster, our fiducial set of 
AGN heating parameters is (if not explicitly stated oth- 
erwise) Ehub = 5 x 10 60 erg, Rhub — 30/i _1 kpc, dbub = 
50/i _1 kpc, and a duty cycle of Atbub = 10 8 yrs, all kept fix 
in time. The thermal energy injected in the form of bubbles 
has been estimated using the simple relations given in Sec- 
tion 2, assuming a~ (5x 10 8 — 3 x 1O 9 )M black hole in the 
cluster centre (depending on the thermal-coupling efficiency 
factor /). For the halos of lower mass, we assumed that the 
thermal content of the bubble is proportional to , in 

analogy to our first scenario for scaling the bubble energy 
in a cosmological setting. This energy scaling is motivated 
by the well established observational relation between the 
black hole mass and the velocity dispersion of the stellar 
component of the bulge (e.g. Tremaine et al., 2002), given 
by 

/ \ 4.02±0.32 

M BH = (1.5 ± O.2)xlO 8 M (^^- rT j , (6) 

and on the hypothesis that the central cluster galaxy scales 
self-similarly with the mass of the cluster itself. Even though 
these assumptions are certainly very restrictive, they provide 
us with a definite model that allows a straightforward inter- 
pretation of the trends with mass, and supply some guid- 
ance for what to expect in full cosmological simulations. We 
also note that the recent numerical work of Di Matteo et al. 
(2003) suggests that, once the Mbh— cr relation is established 
with time, the black hole mass is proportional to of 
the host galaxy. 



3.1 AGN heating of a massive galaxy cluster 

In this section we concentrate on the effects of bubble heat- 
ing on an isolated galaxy cluster of mass 1O 15 /i _1 M , while 
we will discuss the relative importance of AGN heating as 
a function of mass in the next section. In Figures 1 and 2, 
we show maps of the projected mass- weighted temperature 
of the 1O 15 /i _1 M galaxy cluster, focusing on the central 
regions in order to highlight the morphology of the bubbles 
with different injection schemes and at various evolutionary 
stages. Both figures were obtained for simulations with cool- 
ing and star formation, and with AGN feedback of the same 
strength. However, in the left panel of Figure 1, the bubbles 
were placed randomly within a sphere of radius dbub, while 
the remaining three panels illustrate the case of two sym- 



metrical bubbles injected simultaneously, containing half of 
the energy each, and with the injection axis preserved with 
time for different bubble cycles. 

In Figure 2, we show results for simulations with the 
same feedback energy as in Figure 1, but this time the ini- 
tial radius of the bubbles was two times larger and equal 
to 60/i -1 kpc. After being injected, the bubbles rise due to 
buoyancy and start to assume more elongated, "pancake- 
like" shapes, as clearly visible in the left panel of Figure 1. 
They continue to rise until the surrounding gas entropy be- 
comes comparable to their entropy content, at which point 
they settle into an equilibrium and dissolve slowly with time. 
While rising, they push the intracluster gas above them, and 
also entrain some of the cooler central gas, transporting it 
to larger radii. 

A closer comparison of Figures 1 and 2 makes it clear 
that the smaller bubbles with their significantly higher en- 
ergy per particle result in more pronounced mushroom-like 
structures. Nevertheless, they do not shock the surround- 
ing gas which, on the very top of the bubbles, forms cold 
rims. At late evolutionary stages, corresponding roughly to 
a quarter of the Hubble time, (see the right panel of Figure 
2), a characteristic bipolar outflow is visible as a result. 

In Figure 3, we analyse the global gas properties of the 
cluster in terms of radial profiles of density, temperature 
and entropy of the hot gas component, i.e. gas of the cold 
interstellar medium is not included in the plots*. Bubble in- 
jection modifies the inner 100 h~ 1 kpc substantially, reducing 
the density and increasing the temperature profile. Accord- 
ingly, the entropy of the central gas particles is changed as 
well, and an entropy floor is formed. The lower right panel 
of Figure 3 shows the mean gas inflow rate in the central 
30 /i _1 kpc. After a relatively brief period of time, AGN heat- 
ing regulates, in a stable fashion, the flow of gas towards the 
centre, preventing the unrealistically high mass deposition 
rates of a fully developed cooling flow, which can reach up 
to 1200 M yr -1 in the case without bubble heating. Even 
though a repeated injection of bubbles along the same spa- 
tial axis ("jet-like") is somewhat more efficient than a ran- 
dom placement within a sphere, the gas profiles have very 
similar trends in both cases, indicating the robustness of the 



* The 'cold' gas component has been defined here as all gas cooler 
than 1 keV and with a density higher than the star formation 
density threshold. 
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Figure 1. Projected mass- weighted temperature maps of the central regions of an isolated galaxy cluster of mass 1O 15 /i _1 M0. In the 
left panel, bubbles have been introduced with a random placement inside a spherical region, while in the right panel, a "jet-like" injection 
of bubbles is shown where two bubbles are placed opposite of each other, and subsequent generations of bubbles are injected along the 
same spatial axis, i^bubn ^bub an d rfbub i n both cases are the same, and given by 5 x 10 60 erg, 30 h~ 1 kpc and 50 h~ 1 kpc, respectively. 
The maps have been constructed for times of ~ 1.4 Gyr and ~ 0.8 Gyr after the beginning of the runs. 




Figure 2. Time evolution of the isolated 1O 15 /i -1 M0 galaxy cluster with a jet-like AGN heating. The models are the same as in Figure 1, 
only the bubbles have two times bigger radii, namely 60/i -1 kpc. It can be noticed how the morphology in the central cluster region 
changes due to bubble-induced motions, from ~ 1.8 Gyr (left panel) to ~ 3.3 Gyr, when a well-defined bipolar outflow is clearly visible 
(right panel). 



results with respect to these details of the bubble injection 
scheme, at least in situations free of secondary effects due 
to infalling structures and mergers. 

In Figure 4, we show unsharped masked maps of the 
X-ray emissivity of one of our cluster models. The X-ray 
emission has been estimated using the bremsstrahlung ap- 
proximation, 



Lx = 1.2xl0- 24 — Vmg^p,^ 2 [ergs -1 ]. (7) 

LL TTIj tj 

i = 

The unsharp-masking has been performed by subtracting 
from the original projected Lx-map the same map smoothed 
on a 100/i _1 kpc scale. A large number of centrally concen- 
trated ripples are clearly visible in the result. These ripples 
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Figure 3. Radial profiles of gas density (upper left panel), temperature (upper right panel) and entropy (lower left panel) of the isolated 
1O 15 /i -1 M0 halo. Blue lines: run with cooling and star formation. Red lines: AGN feedback mechanism also included (random placement 
of bubbles). The vertical dotted lines denote i?200- Lower right panel: Mean mass inflow rate in the central 30/i -1 kpc as a function of 
time, normalised to the Hubble time. It can be seen that the mass deposition rate onto the central galaxy is substantially reduced with 
AGN heating and stabilised at ~ 15OM0yr -1 after a relatively brief period of time. 



in the X-ray emissivity are sound waves generated by the 
expansion of the bubble after the thermal energy is injected. 
The sound waves travel through the cluster, and if the IGM 
has a residual viscosity they can be dissipated, providing a 
nonlocal heating of the central cluster volume. We note that 
we have explored different scales over which the smoothing 
in the unsharped masked technique is performed, obtain- 
ing the sharpest and most prominent features for smoothing 
scales corresponding to approximately 100 /i _1 kpc, which by 
order of magnitude agrees with the dimension of the ripples 
themselves. 

We find that the ripples reach distances of ~ 800 h~ 1 kpc 
after 1 Gyr, translating to a velocity of order ~ 10 3 kms -1 , 
which matches the expected sound speed in the ICM of this 
cluster. At larger radii, the sound waves are not detectable 
any more. Note that we also expect that their velocity drops 



strongly in the outskirts of the cluster, where the tempera- 
ture and hence the sound speed decline. 

Upon closer inspection, it can be seen that the ripples 
are actually slightly offset from the cluster centre, with their 
midpoint directly matching the initial coordinates of the in- 
jected bubble. Moreover, the ripples progressively lose their 
intensity at larger radii, both due to a 1/r 2 dilution of their 
intensity, and to a lesser extent, due to a damping caused 
by the residual viscosity of our SPH scheme. Note that some 
level of numerical viscosity is instrinsic to all SPH schemes, 
even though we are modelling an ideal gas. Quantifying 
the exact magnitude of the resulting effective viscosity is 
not trivial, also because it depends on the spatial resolu- 
tion achieved in the simulations. However, recent observa- 
tions of optical Ha- filaments (Fabian et al., 2003b) suggest 
that the gas in the central region of the cool-core cluster 
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Figure 4. Unsharped masked map of the X-ray luminosity in the 
central region of the 10 15 h~ 1 M.Q isolated halo, at time ~ 2.2 Gyr. 
The unsharped masking has been performed by subtracting a 
smoothed map from the original projected X-ray emissivity map, 
with the smoothing scale set to 100 /i -1 kpc. It can be seen that 
the AGN bubble heating generates a number of sound waves, 
which could gradually release their energy to the ICM if they are 
viscously damped on their way to the cluster outskirts. 



Perseus might be quite viscous, rather then turbulent (but 
see Ensslin & Vogt, 2005, for estimates of gas turbulence on 
smaller scales). If ICM viscosity is relevant, it would imply a 
high rate of dissipation of the energy contained in the sound 
waves at small radii. This physical viscosity could be easily 
higher than the numerical viscosity we have in our simu- 
lations. Naturally, it is then desirable to treat the dissipa- 
tion process accurately, which requires an SPH discretization 
of the Navier- Stokes equation combined with an assumed 
level of physical Spitzer- viscosity. Recently, first mesh-based 
studies of isolated clusters with viscosity and bubble heat- 
ing have appeared (e.g. Ruszkowski et al., 2004; Reynolds 
et al., 2005; Briiggen et al., 2005). We plan to investigate 
this theoretical issue in a forthcoming study. 

In Figure 5, we show the locus of selected particles in 
the log S — log R plane, at a time equal to one quarter of the 
Hubble time (which marks the end of most of our isolated 
simulations). We selected only particles that at least once 
belonged to one of the injected bubbles. For easier compari- 
son, the mean radial entropy profile of the AGN heated clus- 
ter (the same as in Figure 3) is plotted as a red dot-dashed 
line. The dots shown for each individual bubble particle have 
been colour- coded according to their relative temperature. 
The particles of a recently injected bubble have the high- 
est temperature values. As expected, their entropy values 
lie substantially above the average entropy of the cluster for 
the same range of radii, implying that the bubble will rise 
due to buoyancy. We thus expect these bubble particles to 



reach larger radii and to lose some of their thermal energy 
content during the rise, and this expectation is borne out by 
the cooler particles from older bubbles. 

The radius at which the bubble particles attain the 
mean cluster entropy level is set by their initial density and 
thermal content after injection, by their capacity to shed 
some of the energy to the surrounding ICM, by their radia- 
tive cooling efficiency, and by the amount of mixing. Given 
that the different bubbles at various epochs have very simi- 
lar initial temperature and mass, Figure 5 implies that the 
bubbles reduce their temperature almost by one order of 
magnitude, from the injection instant to the final equilib- 
rium position. If we sum up the total energy injected over 
the entire simulated time, and assume that it gets thermal- 
ized over the whole cluster, we obtain that the gas particle 
temperature is increased by ~ 0.2 keV, which roughly cor- 
responds to the bubble temperature decrement mentioned 
above when the mass fraction of the bubbles is taken into 
account. Thus, it appears that the radiative cooling is not 
severe inside the bubbles, even though in the cluster as a 
whole it approximately balances the AGN feedback mecha- 
nism. 

We extended our investigation by considering 'jet-like' 
injection of bubbles and also a scenario in which the bubbles 
are inflated in a continuous fashion over some time interval 
tinj • We tried values from ti n j = 5 x 10 7 yrs to ti n j = 5 x 10 8 yrs, 
which is significantly longer than the sound crossing time 
over the scale Rbub of the bubble. The maximum radius 
reached by the bubbles is essentially invariant in all of these 
cases, yielding ~ 250 — 300 /i -1 kpc at the final simulated 
epoch. Also, the heating efficiency of buoyant bubbles re- 
mains very similar, although the bubbles are somewhat more 
energetic in the continuous injection scheme, presumably be- 
cause cooling losses are reduced here due to the expansion of 
the bubble before the bulk of the energy us released. This is 
directly reflected in an even lower mass deposition rate onto 
the central object, which always occurs in a stable fashion 
with time, where the gas cooling inflow is balanced by the 
AGN heating rate. 

It is important to point out that the entropy content 
of the bubbles and the maximum distance they can reach 
from the cluster centre depend upon the equation of state 
assumed for the gas belonging to the bubbles. In all our 
models bubbles have been simulated assuming the equation 
of state of an ideal gas. However, radio observations indi- 
cate that AGN-driven bubbles contain relativistic particles 
which possibly dominate over the thermal pressure compo- 
nent, implying a softer equation of state. Moreover, the en- 
ergy contrast of the individual bubbles is not very high in our 
approach, resulting in a gentle ICM heating, without pres- 
ence of significant shocks. In fact, most of the observations of 
AGN-heated clusters point out that strong bubble-induced 
shocks appear to be absent, although recently a few clusters 
with moderate shocks in connection with AGN activity have 
been discovered (Fabian et al., 2003a; Nulsen et al., 2005a; 
McNamara et al., 2005). Therefore, due to the assumptions 
of our model, the maximum possible distance reached by 
the buoyant bubbles in the cluster atmosphere may be un- 
derestimated with respect to the case where the relativistic 
particle component is modelled as well. 

An important question is whether the bubbles are ca- 
pable of raising cold gas from the cluster centre and mixing 
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Figure 5. Mean radial entropy profile of the 10 15 h~ 1 M.Q isolated 
halo. The result for the simulation with AGN feedback is given 
by the red dot-dashed line. The dots show the positions of the 
bubble particles, and they are colour-coded according to their 
temperature. 



it with higher entropy gas at larger radii. Note that the X- 
ray observations of central metal abundance gradients put 
a constraint on the amount of gas mixing (Bohringer et al., 
2004) in the centre. In order to address this issue, we anal- 
ysed the gas metallicity distribution in the clusters. With- 
out AGN feedback mechanism, all metals produced in the 
10 15 h~ 1 Mo isolated halo are enclosed in the star forming 
region and are confined to the very centre, where the density 
is sufficiently high to allow star formation. The bubble heat- 
ing instead produces both a reduction of the star formation 
rate by heating some of the gas that would otherwise end-up 
in the cold interstellar medium, and a spreading of metals 
away from the cluster centre. Moreover, in our simple model, 
the metals produced by the central stars are partly entrained 
and transported along with the bubbles to larger radii. The 
spatial distribution of the stars themselves is unaffected by 
the bubbles, however. It is important to note that the metal 
mixing in our model due to the bubbles represents a lower 
limit on the induced additional mixing, because fluid dy- 
namical processes that produce small-scale mixing tend to 
be underresolved in cosmological simulations. 



3.2 Efficiency of bubble heating in halos of 
different mass 

The radiative cooling times of halos depend on their mass, 
both because of their different virial temperatures, and be- 
cause of the temperature dependence of the cooling rate. 
Given that the masses of supermassive black holes, and 
hence our assumed bubble feedback, have a mass depen- 
dence as well, we expect a complex interplay between the 



different heating and cooling processes, and as a result a 
mass-dependent efficiency of the feedback. This non-linear 
dynamics can be best studied with detailed numerical anal- 
ysis. To this end, we simulated isolated halos for a range 
of masses, starting from 10 12 /i _1 M© and reaching up to 
10 15 h~ 1 M ( 7 ) , with the characteristics listed in Table 1. The 
case of the most massive halo has been discussed in some de- 
tail in the previous section, so that we can restrict ourselves 
here to highlight the differences that occur when the mass 
of the systems is lowered. The study of smaller halos gives 
also direct insight into the question of how bubble feedback 
may affect the hierarchical assembly of present-day massive 
clusters. 

In our numerical simulations, the coupled dynamics re- 
sulting from cooling, subsequent star formation and a given 
heating mechanism is quite complex. The introduction of 
a certain amount of heating can in special situations even 
trigger an increase of the net cooling rate. For example, let 
us consider star-forming gas with a density higher than the 
density threshold set for star formation. If this cold gas com- 
ponent receives an amount of thermal energy from bubble 
heating that is insufficient to bring it back into the hot 
phase (where most of the intracluster gas resides), then a 
counterintuitive process may occur. In this case, the ther- 
mal energy injection prevents the cold gas component from 
forming stars, but the local gas density will remain compar- 
atively high, which in turn stimulates even larger radiative 
cooling losses. Thus, such a gentle heating, especially in low 
mass systems where radiative cooling is more pronounced, 
can in extreme cases even stimulate an increase of the cold 
gas component. 

Analysing diagnostic phase-space diagrams like the 
log p — logT plane, one can notice that for the 10 15 /i _1 M© 
halo the relative quantity of central cool gas is low, and 
it is promptly heated once bubble injection is switched on. 
Moreover, the energy per bubble particle is sufficiently high 
to push the gas into the "hot-phase" (upwards and to the left 
in the diagnostic diagram, as illustrated for the 10 14 /i _1 M© 
halo in Figure 6) and the star formation at late simulated 
epochs is completely quenched. 

Examining the 10 14 h~ 1 M.Q isolated halo, we find that 
bubbles are still very efficient in reducing the star formation 
rate, e.g. after the time £h/4, only 14% of stars are formed 
with respect to the run without AGN heating, but the total 
amount of cold gas, both in the central regions and out to 
i?2oo, remains very similar. The diagnostic diagram for this 
cluster is shown in Figure 6. The small orange dots are the 
gas particles for the run without AGN feedback, and they 
maintain practically the same position even when bubble 
heating is included. The big blue dots are those gas parti- 
cles that belong to the 'cold phase', here defined as having 
temperatures less than 1 keV and densities higher than the 
density threshold set for star formation (as indicated by the 
vertical dotted line). With AGN heating included, the red 
dots denote the bubble particles, while the green star sym- 
bols give the locations of cold phase gas particles. Two differ- 
ent features due to the presence of bubble feedback are read- 
ily apparent. First, the cold gas fraction for the intermediate 
temperature range, from 5 x 10 5 K to 10 7 K, is substantially 
reduced in the case with feedback, because for most of these 
particles it is possible to transport them back to the hot 
phase. In contrast, along the line in the lower-right part of 
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Figure 6. Phase-space diagram of gas temperature versus gas 
density for the 10 14 h~ 1 M.Q galaxy cluster. Small orange dots 
are the particles outside the cold star-forming region, while big 
blue dots denote the gas particles in the run with no feedback, at 
high overdensities and with temperatures below 1 keV. Green star 
symbols are the particles satisfying the same criteria, but when 
the AGN heating is included. Finally, the position of the bubble 
particles is given by red dots. 



the diagram, there are more particles when bubble heating 
is included. This can be explained by the fact that this line 
is determined by the multiphase structure of the ISM, given 
here in terms of an effective mass-weighted mean tempera- 
ture, which is a combination of the temperature of cold gas 
clouds and the one of the hot ISM component (see Springel 
& Hernquist, 2003). Hence, while the feedback mechanism 
reduces the number of stars formed, it produces a higher 
amount of interstellar gas with very low temperatures. Note 
that a large fraction of these cold gas particles have been a 
part of a bubble at some earlier epoch. 

In the lower mass system of mass 1O 13 /i _1 M0, the fi- 
nal number of stars is also reduced, this time by only 13%, 
and, as well as in the previous case, the cold gas fraction 
is essentially unchanged. In the most extreme case of the 
10 12 /i _1 M© halo, the bubble heating of the ICM is radi- 
ated away on a very short timescale without producing any 
substantial modification. 

We considered different injection assumptions in order 
to test whether this effect can be alleviated and the bubble 
heating efficiency can be increased. Focusing on the most dif- 
ficult case of the 10 12 /i _1 M© halo, we explored a range of in- 
jection energies, from 5 x 10 56 erg per bubble, to 3 x 10 57 erg. 
Also, we injected bubbles according to different schemes: in a 
spatially correlated 'jet-like' fashion, or by inflating the bub- 
bles gradually for ti n j = 5 x 10 7 yrs, i.e. by releasing the en- 
ergy slowly instead of instantaneously. Finally, we also tried 
a model where we artificially prevented bubble particles from 
cooling for 10 7 — 10 8 yrs, motivated by some observational ev- 



idences that the bubbles contain a non-thermal component, 
which would then be able to maintain its pressure longer. 
From these experiments we conclude that relating the bub- 
ble energetics to the underlying dark matter potential by 
invoking an assumed relation with its black hole mass, does 
not in general provide stable solutions for an efficient elimi- 
nation of cooling flows in low mass systems. Unrealistically 
high bubble energies are required in order to offset the cool- 
ing flow, and also their thermal content has to be fine-tuned 
to a restricted range of values, otherwise AGN heating may 
easily become so strong that bubbles blow out substantial 
amounts of mass from the halo or cluster potential well. 

This clearly indicates an important deficiency of bub- 
ble models like the one studied here. Since there is no in- 
ternal trigger for bubble activity, a self-regulation loop is 
missing, but this will be required for stability in more gen- 
eral situations. Therefore, a more detailed physical scenario 
for the triggering of bubble activity is needed which couples 
the local physics of the cooling flow with the activity and 
energetics of AGN feedback. We will discuss a number of 
possibilities for this in Section 5. 



3.3 Observational X— ray features of simulated 
bubbles 

In the last few years, a growing number of observations per- 
formed with the Chandra X-ray telescope have found the 
presence of so called X-ray cavities in central galaxy clus- 
ter regions (e.g. McNamara et al., 2000; Fabian et al., 2000; 
Sanders & Fabian, 2002; Mazzotta et al., 2002; Birzan et al., 
2004). These observations support a scenario where relaxed 
cooling-flow/cooling-core clusters are heated due to the pres- 
ence of a supermassive black hole at their centres. Hence, it 
is interesting to see whether the simulated bubbles produced 
by our model have morphologies comparable to the observed 
ones. 

In order to perform this comparison accurately, it is 
necessary to produce artificial emissivity maps which are as 
similar as possible to realistic observations based on a finite 
exposure time on an X-ray telescope. Recently, a similar 
analysis has also been performed by Briiggen et al. (2005). 
To this end, we processed selected simulation outputs with 
the X-MAS software package. A final result of this code is a 
photon event file quite similar to the one an observer would 
acquire with the Chandra telescope in ACIS-S3 mode. The 
instrument background has been included in our images by 
taking into account an appropriate blank-sky background 
file (Markevitch et al., 2001). A detailed description of the 
X-MAS simulator can be found elsewhere (e.g. Gardini et al., 
2004); here we limit the description to the subsequent anal- 
ysis steps performed and the significance of the maps ob- 
tained. 

We have generated event files for different exposure 
times, ranging from 10 ks to 1 Ms, both for the runs with 
and without AGN feedback. We then selected a number of 
different energy bands, performed a Gaussian smoothing on 
a range of scales, and finally produced unsharp masked im- 
ages to search for evidence of systematic departures of the 
flux from the mean. In Figure 7, we show photon images of 
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the central region of the 10 15 /i _1 M© isolated halo^, both in 
a case with and without additional AGN heating. The phys- 
ical scale of the maps corresponds to ~ 670 kpc (2048pix), 
the energy band has been chosen to be AE = [0.3, 1.5] keV, 
and the maps have been smoothed by summing the pixel 
fluxes in bins of 4 pixels. For this AE, the instrument back- 
ground is minimised and the features due to the presence of 
the bubbles are more evident. 

The first panel of Figure 7 shows a photon image of 
the AGN-heated cluster after 100 ks of exposure time, be- 
fore applying any smoothing. The rest of the plots have been 
created by Gaussian smoothing them first on a small scale 
(3pix), then re-smoothing the obtained image on a bigger 
scale (15pix), and finally unsharp masking the two smoothed 
images (i.e. subtracting off the 15pix smoothed version). The 
smoothing scales have been selected to maximise flux depar- 
tures from the mean. The second and the third panel illus- 
trate how the bubbles introduce emissivity irregularities for 
two different exposure times, 100 ks and 1 Ms, respectively. 
Finally, the fourth panel presents the cluster photon image 
after an exposure time of 100 ks and with no AGN feedback. 

It is clear that the bubbles generate characteristic fluc- 
tuations in the photon counts, both creating bright features 
and X-ray depressions. The typical dimension of these irreg- 
ularities is ~ 50 kpc, very similar to the size of the bubbles 
themselves. The hot spots can be associated with the most 
recent bubble events, containing particles still significantly 
hotter than the surrounding ICM (as can be also seen from 
Figure 5), whereas the depressions in photon counts can be 
explained with previous bubble episodes. These peculiari- 
ties in emissivity are completely absent in the galaxy cluster 
without AGN feedback, indicating that they are real features 
and not artifacts produced by counting statistics or our anal- 
ysis. The only feature that is present in the fourth panel of 
Figure 7 is the central excess due to the prominent cooling 
flow of rsj 100 kpc size in diameter. Based on these results we 
conclude that in a relaxed galaxy cluster, departures from 
the mean flux stemming from bubbles with characteristics 
as given by our model can be detected, provided the expo- 
sure times are long enough, and provided that other sources 
of photon count deviations are absent or negligible. 



4 EFFECT OF AGN BUBBLE HEATING IN 
COSMOLOGICAL SIMULATIONS 

4.1 Simulation characteristics 

As a next step in our analysis, we consider the importance 
of AGN feedback in full cosmological simulations of cluster 
formation. To this end, we selected a number of galaxy clus- 
ters with a wide range of masses from a parent dark matter 
simulation, and resimulated them with higher resolution in- 
cluding gas dynamics. Our hydrodynamical simulations ac- 
count for cooling and star formation, and additional AGN 
heating as well. In a subset of our runs, we also included 

t We decided to perform this analysis on an isolated cluster in or- 
der to minimise other features in the X-ray emissivity that would 
have been imprinted by possible substructures or merger events. 
Further discussion of this issue in the cosmological framework is 
given in Section 4.5. 



galactic winds powered by star formation, as implemented 
by Springel & Hernquist (2003). Models with winds pro- 
vide a better description of some galaxy cluster properties, 
e.g. the distribution of metals, but in order to be able to 
cleanly identify the effects of bubble heating, we focus most 
of our analysis on simulations without winds. Where appro- 
priate, we will however briefly discuss any changes of our 
results when winds are also included. 

Our primary series of simulations consist of resimula- 
tions of a cluster extracted from the GIF ACDM simula- 
tion (Kauffmann et al., 1999). We selected the second most 
massive galaxy cluster from this simulation and constructed 
higher resolution initial conditions for it using the "Zoomed 
Initial Conditions" technique (Tormen et al., 1997). We car- 
ried out two runs with different resolution in order to test nu- 
merical convergence of our results. These clusters are equiv- 
alent to the ones used by Springel et al. (2001a), but with 
gas (from now on we will refer to these runs as SI and S2, 
respectively). The simulations have been evolved from an 
initial redshift of Zi n i = 30 for SI, and Zi n i = 50 for S2, 
producing 25 outputs uniformly spaced in the logarithm of 
the expansion factor. Additionally, we selected two other 
clusters with final virial masses substantially smaller (g676) 
and bigger (gl) than the Sl/S2-cluster. These clusters have 
been extracted form a cosmological ACDM simulation of 
box- size 479/i _1 Mpc (Yoshida et al., 2001; Jenkins et al., 
2001), and again have been resimulated with the ZIC tech- 
nique at higher resolution (Dolag, 2004). 

The simulation of the gl galaxy cluster includes sev- 
eral other smaller systems in the high-resolution region 
which we also included in our analysis. Tables 2 and 3 pro- 
vide a summary of the main properties of our set of simu- 
lated galaxy clusters. In all runs, the cosmological param- 
eters were that of a flat concordance ACDM model, with 
Qm = 0.3, Qa = 0.7, Qb = 0.04, a normalisation of the 
power spectrum given by as = 0.9, and a Hubble constant 
at the present epoch of H — 70 km s _1 Mpc -1 . 

Unlike simulations of isolated clusters, cosmological 
simulations require a special algorithmic method for placing 
bubbles, since the position of the cluster centre and proper- 
ties like virial mass are not known a priori, and change with 
time. To address this problem, we run for every AGN duty 
cycle a fast parallel FOF group finder on the fly as a part of 
the simulation code, obtaining a list of all halos with their 
basic properties. We then adopt two different schemes for 
injecting bubbles. We either consider only the most massive 
halo found in the high-resolution zone, which can be identi- 
fied with the most massive progenitor of the final cluster, or 
we introduce AGN-driven bubbles in all halos above a given 
fixed mass threshold value. The injection of bubbles in all 
large halos is motivated by the observational indications that 
probably most if not all of the spheroidal galaxies harbour 
a supermassive black hole at their centres. Note that the 
larger number of bubbles in this second scenario can also 
cause additional effects during merger events, where bubble 
material can be torn apart and mixed into outer regions of 
the cluster. 
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Figure 7. Artificial photon images for the 10 15 h~ 1 MQ isolated galaxy cluster obtained with the X-MAS software package. In all the 
panels, photons have been selected in the energy band AE = [0.3, 1.5] keV, maps were binned in 4pix/bin, and the physical scale of the 
maps is ~ 670 kpc. The cluster emissivity map before applying any smoothing is illustrated in the first panel, while all the other panels 
have been obtained by unsharp masking previously smoothed images, as explained in the text. The second and the third panels show 
qualitatively very similar features. They are for the run with AGN feedback and differ only in the exposure time. Finally, the fourth 
panel shows how the same cluster appears when the bubble heating is absent. 



Simulation 



A^HR 



At 



gas 



m D M [h 1 M Q 



[h 1 M ] Zstart Zend € [h X kpc] 



51 450088 450088 

52 1999978 1999978 
g676 314518 314518 

gl 4937886 4937886 



5.96 x 10 9 
1.18 x 10 9 
1.13 x 10 9 
1.13 x 10 9 



0.92 x 10 9 
0.18 x 10 9 
0.17 x 10 9 
0.17 x 10 9 



30 
50 
60 
60 



14.5 
8.5 
5.0 
5.0 



Table 2. Numerical parameters of the cosmological galaxy cluster simulations used in this study. The values listed from the second to 
the fifth column refer to the number and to the mass of high resolution dark matter particles and of gas particles. Note that the actual 
values of A^ gas and m ga s vary in time due to star formation. The last three columns give the initial and final redshifts of the runs, and 
the gravitational softening length e. 



4.2 Global gas properties of simulated galaxy 
clusters 

Before analysing the properties of simulated galaxy clusters 
with and without AGN bubble heating, we briefly discuss is- 
sues of numerical convergence. For this purpose we consider 
the SI and S2 runs, and compare their spherically averaged 
radial profiles at two epochs, namely at z — 3 and z — 0*. 
The dark matter and stellar density profiles of the SI and S2 
galaxy clusters are in excellent agreement at both epochs, 
as well as the gas density profiles, with the residual differ- 
ences at early times being consistent with what is expected 
from the increased noise. Both the emission- weighted and 
the mass-weighted temperature profiles do not noticeably 
differ at low redshifts, while there is a hint a of slightly higher 
gas temperature for the S2 cluster at early times. Thus, we 
conclude that for radii larger than the gravitational soften- 
ing length the properties of our simulated galaxy clusters 
are numerically robust and have converged quite well. 

In Figure 8, we compare the gas entropy profiles with 
and without bubble heating at three different epochs, z = 
1.64, z = 0.44 and z = 0. For this comparison, we use both 
of our AGN heating models, the one based on the Magorrian 
relationship and the "BHAR model" . The entropy has been 
estimated by calculating the ratio of the emission-weighted 
temperature to the gas density elevated to the 2/3 power, 



$ These two epochs delimit the time interval during which the 
bubble heating is active, and hence the period of time where our 
analysis is performed. 



where the temperature is measured in Kelvin and the gas 
density is given in /i 2 M0kpc -3 . We selected only the hot 
gas component to compute these profiles, i.e. we avoided the 
cold, star- forming gas by imposing a cut in density and ion- 
isation level. With this choice, the gas profiles are smoother 
because they have no contributions from cool substructures 
at various radii. Nonetheless, it is also important to inves- 
tigate the fate of the cold gas in the central cluster region, 
an issue we will address separately in Section 4.3. The blue 
continuous lines are for the run without AGN feedback, the 
red dot-dashed lines correspond to the "Magorrian model", 
while the green dashed lines are for the "BHAR model". 
The vertical dotted lines denote the softening length and 
the virial radius at the different epochs, respectively. When 
^bub is computed from the "BHAR model", the effect of 
bubbles is less prominent at low redshifts than in our other 
AGN heating scenario. However, at early times the situation 
is opposite, as expected. Here the "BHAR model" heats the 
ICM gas more prominently, right from the initial injection 
epoch (z — 3) until z ~ 0.4. It is interesting to note that 
at z ~ 0.4 the bubble energy content is already much lower 
than in the "Magorrian model" , indicating that the efficient 
heating at early times has a prolonged effect on the ther- 
modynamic state of the galaxy cluster. We find that the 
"Magorrian model" starts to affect the gas entropy from 
z « 1.6 in a noticeable way and it becomes very important 
at late times, where it suppresses the cooling flow completely 
at z — 0. As a further comparison, we show in Figure 9 ra- 
dial profiles of gas density, temperature, X-ray luminosity, 
and local cooling time for the SI galaxy cluster at z — 0. 
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Cluster #200 [/i _1 kpc] M 2 oo[^" 1 M ] T mw [K] T ew [K] Lxlergs" 1 ] 



SI 


2427 


9.98 x 10 14 


5.2 x 10 7 


8.7 x 10 7 


9.8 x 10 44 


S2 


2466 


1.05 x 10 15 


5.1 x 10 7 


8.8 x 10 7 


9.9 x 10 44 


g676 


1176 


1.13 x 10 14 


1.4 x 10 7 


2.6 x 10 7 


1.6 x 10 43 


gl_a 


2857 


1.63 x 10 15 


7.3 x 10 7 


1.3 x 10 8 


1.0 x 10 45 


gl-b 


1914 


4.89 x 10 14 


3.1 x 10 7 


4.1 x 10 7 


1.2 x 10 44 


gl_c 


1448 


2.12 x 10 14 


1.5 x 10 7 


2.5 x 10 7 


3.0 x 10 43 


gl_d 


1258 


1.39 x 10 14 


1.6 x 10 7 


1.9 x 10 7 


1.8 x 10 43 


gl_e 


1085 


8.92 x 10 13 


1.1 x 10 7 


1.6 x 10 7 


7.0 x 10 42 



Table 3. Physical properties of our sample of simulated galaxy clusters at z = and at 200p c . For different galaxy clusters, labeled in 
the first column, cluster radius, total mass, mass- and emission-weighted gas temperature and X-ray luminosity are listed, respectively. 
Note that the values refer to the simulations with cooling and star formation, without bubble heating included. 




Figure 8. Radial profiles of gas entropy of the SI galaxy cluster simulation. The blue continuous lines are for the run with cooling and 
star formation only, the red dot-dashed lines refer to the case when AGN heating based on our "Magorrian-like scheme" is included as 
well, while the green dashed lines are for the BHAR-based model. The dotted vertical lines denote the gravitational softening and the 
virial radius at the given redshift; the latter is indicated in the upper-left corner. The profiles do not extent down to vanishingly small 
radii because they have been calculated exclusively from the hot gas component (basically the gas above 1 keV), excluding the cold dense 
gas in the centre. Note that both the spatial and the entropy scale vary between the three different panels. 



AGN heating alters the gas properties out to a radius of 
« 300/i _1 kpc, reducing the central gas density and increas- 
ing its temperature. The X-ray emissivity, being more sen- 
sitive to the gas density, is substantially lower when bubble 
feedback is active. In the lower-right panel of Figure 9, we 
show the cooling time of all ICM gas, estimated isobarically 
as (Sarazin, 1988) 

^' = 8 - 5xl0 iiF^) _1 (i^K) 1/2 ^' < 8 > 

where n p is the number density of hydrogen. When AGN 
heating is not included, the cooling radius, i.e. the radius 
where t coo i = ^h, lies at ~ 60/i _1 kpc, while it gets reduced 
to - 25/i _1 kpc in our "BHAR model". However, the cool- 
ing radius vanishes for the "Magorrian model", where the 
bubbles injection heats the gas above 1 keV. 

Even though the spatial extent of the bubble particles 
reaches out to the virial radius of the cluster, they are not 
capable of heating the gas in outer regions, simply because 
their entropy content becomes comparable to the entropy of 
the surrounding ICM gas at intermediate radii. This find- 
ing is analogous to the previously discussed case of isolated 



halos, but the dynamical evolution of the cluster with its 
associated merger processes makes spreading of bubble ma- 
terial towards the outskirts more efficient. 

In Figure 10, we show emission- weighted temperature 
maps of the g676 galaxy cluster at four different epochs, in 
order to more closely discuss the spatial distribution of bub- 
bles during merger events. The over-plotted dots represent 
the particles that at least once belonged to a bubble, and 
they are colour- coded according to their temperature, the 
darkest ones are particles with T > 10 8 K, while the lightest 
have T < 10 4 K. For this analysis, we have introduced AGN 
feedback in all halos above 5 x 10 10 /i _1 M© and we scaled 
the energy content of the bubbles with M^^z) of the host 
galaxy cluster. 

In the first panel at redshift z — 0.86, there is a smaller 
halo on the lower right corner which enters the most massive 
cluster progenitor at that epoch (which roughly has three 
times larger mass), which is located at the centre of the 
panel. Both the massive halo and the smaller one are AGN 
heated, but the bubbles are less energetic for the infalling 
halo due to our assumed mass dependence. Moreover, it can 
be noticed that the bubbles occupying the central regions 
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Figure 9. Radial profiles at z = of gas density (upper-left panel), emission- weighted temperature (upper-right panel), and X-ray 
luminosity (lower-left panel) estimated with the bremsstrahlung approximation given by eq. (7). The lower-right panel shows the cooling 
time of all gas particles as a function of radius, computed using eq. (8). The continuous horizontal line indicates the Hubble time at z = 0. 
The blue continuous lines are for the case without AGN feedback, the red dot-dashed lines are for the model where i^bub °c M^qq (2)5 
while the green dashed lines are for the scenario where the bubble energy depends on the BHAR given by eq. (3). 



are hotter, both because they are more recent and thus have 
had less time to lose their energy content, and also because 
at later times bubbles are intrinsically more energetic in our 
"Magorrian scenario". The second panel of Figure 10 (at 
z — 0.62) illustrates what happens to the bubble distribu- 
tion when the smaller halo is crossing the central region of 
the massive cluster. The bubbles are literally pushed out of 
the way, upwards and to the right, and they are also heated. 
The next panel (at z — 0.48) shows that the bubble dis- 
tribution is still quite asymmetric, but at the same time, 
bubbles have spread efficiently into outer regions and also 
have cooled. Finally, the last panel shows how the cluster 



appears at z — 0.20, where it starts to be fairly relaxed with 
a quite symmetric distribution of bubbles. 

4.3 Stellar properties of galaxy clusters 

In this section, we analyse the effects of AGN heating on the 
properties of the stellar components of galaxy clusters. We 
concentrate on the properties of the central cluster galaxy, 
which is the one affected most by the bubble heating. From 
the initial injection epoch (z = 3) until z = 0, we compute 
for this purpose the stellar and gaseous mass, star formation 
rates, stellar ages, and colours of the cD galaxy that sits in 
the main progenitor of our final halo. 
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Figure 10. Emission- weighted temperature maps of the g676 galaxy cluster simulation during a major merger event at z = 0.86, 0.62, 
0.48 and 0.20, respectively. Over-plotted dots represent gas particles that have been at least once part of a bubble, and they are colour- 
coded according to their temperature, the darkest ones being the hottest. It can be noticed that both the spatial distribution of bubbles 
and their energy content are drastically influenced by the passage of the smaller halo through the central region of the massive cluster. 



In Figure 11, we show a histogram of the formation 
times of stars belonging to the cD galaxy at z = of our SI 
galaxy cluster. The histograms have been computed by bin- 
ning the expansion factors that correspond to each stellar 
formation time, and they have been normalised to the max- 
imum bin. The white histogram is for the run without AGN 
heating, the grey coloured one for the "Magorrian model", 
while the hatched histogram gives the result for "BHAR 
model". When AGN heating is absent, the histogram of stel- 
lar formation redshifts clearly shows an extended tail at low 
redshifts; together with the considerable SFR of order of 
100M©yr -1 , this implies that stars are formed until z — 
in situ. Nevertheless, there is a possibility that some of these 
stars have been formed elsewhere, e.g. in merging substruc- 



tures, and that they only ended up later in the cD galaxy by 
merging. For redshifts less than 0.3 this is certainly not an 
important mechanism, because in our "Magorrian model", 
the central SFR is completely suppressed for z < 0.3 and 
at the same time the z s f histogram is truncated. Thus, the 
difference of ~ 6 x 10 11 h~ 1 Mo in the stellar mass of the fi- 
nal cD galaxy in these two cases gives an indication on how 
many stars have formed in the cD galaxy from z = 0.3 un- 
til today, when AGN feedback is not present. At the same 
time, also the mass of the cold gas (below 1 keV) is reduced 
in the "Magorrian model", from ~2x 10 10 /i _1 M© to zero. 
Moreover, the suppression of star formation at late times 
has an immediate effect on colours of the cD galaxy, which 
becomes redder. To estimate the colours we used Bruzual 
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& Chariot's stellar population synthesis models (Bruzual 
& Chariot, 2003), computing rest-frame magnitudes in the 
SDSS bands, assuming Solar metallicity and a Chabrier ini- 
tial mass function. The i^-band magnitude is changed from 
—23.8 to —22.8 and the u — r colour is increased from 2.0 to 
2.6. 

In contrast, the star formation of the "BHAR model" 
proceeds in a quite different manner, as can been seen 
from the hatched histogram, which is substantially lower 
for 0.4 < z < 2.0, but quite similar to the case without 
AGN feedback at very low redshifts, z < 0.3. The first fea- 
ture arises due to two different processes happening at the 
same time. The second peak present in the histograms is 
due to a major merger event that happens roughly at z = 1. 
One consequence of this merger event is a central burst of 
star formation, which happens to be absent in the "BHAR 
model" due to its very efficient bubble heating at this epoch. 
Nonetheless, there are still some stars becoming part of the 
cD galaxy at z = that have formation times in the corre- 
sponding time interval. Tracing back these stars in time and 
considering that the SFR of the cD galaxy in the "BHAR 
model" during this epoch is practically zero, we see that 
these stars have been produced in other small galaxies and 
were indeed accreted onto the central cD galaxy at later 
times. Assuming that a similar amount of stars accreted 
onto the cD galaxy also in the case without AGN feedback, 
a total mass of ~ 10 12 /i _1 M© in stars has formed in situ 
between z — 2 and z — 0.4, with a considerable part created 
as part of the central starburst induced by the major merger 
event. 

The tail at low redshifts present in the hatched his- 
togram can be explained by the reduced energy content of 
the bubbles, which are not efficient any more in suppress- 
ing the cooling flow, and consequently the star formation in 
the central cD galaxy. The total SFR within R200 follows 
a similar trend as the SFR of the cD galaxy, with compa- 
rable systematic differences for the different runs, implying 
that the bubble heating mainly affects the central stellar 
properties and not the residual star formation in the cluster 
volume. 

4.4 The metallicity distribution in the simulated 
clusters 

In this section, we analyse the metal distribution in our sim- 
ulated galaxy clusters. It is a well known problem that nu- 
merical simulations which include cooling and star formation 
processes in general fail to reproduce the observed shallow 
metallicity gradients, especially so if efficient feedback mech- 
anisms that help spreading the metals are absent. In partic- 
ular, metals produced by stars remain locked in the dense, 
star forming regions, even though supernovae feedback is in- 
cluded which regulates the star formation process itself. As 
a result, the metallicity distribution remains lumpy and ex- 
hibits a rather step gradient. Moreover, most of the metals 
are produced in the central cD galaxy due to its excessive 
star formation, which is a manifestation of the cooling flow 
problem. 

This motivates the search for physical feedback pro- 
cesses that can spread and mix metals more efficiently, act- 
ing both in the central region and on the scale of the whole 
galaxy cluster. While the galactic wind model suggested by 
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Figure 11. Distribution of formation redshifts of the stars be- 
longing to the cD galaxy of the SI galaxy cluster at z = 0. The 
white histogram corresponds to the case with cooling and star for- 
mation only; here it can be seen that stars continue to be formed 
until z = 0. The grey coloured histogram shows how the stellar 
formation times change when our "Magorrian model" of bubble 
heating is switched on, which essentially suppresses any central 
star formation for z < 0.25 completely. Finally, the result for the 
"BHAR model" of AGN feedback is illustrated with the hatched 
histogram, which shows a considerable suppression of central star 
formation at intermediate redshifts, 0.3 < z < 2.0. 



Springel & Hernquist (2003) helps in reducing this discrep- 
ancy and also diminishes the total SFR over cosmological 
time, the model fails to qualitatively change the gas and 
stellar properties of galaxy clusters discussed in Sections 4.2 
and 4.3. Especially at late times, mixing of metals due to 
winds from the cD galaxy proves inefficient; here the cluster 
potential well is simply to deep and the ram pressure of the 
ICM too high to allow winds to travel far. Bubble heating 
may fare considerably better in this respect. In the following 
we therefore analyse the effect of AGN on the metal distri- 
bution in our simulations, and we also compare simulations 
with or without galactic winds of velocity ~ 480kms _1 . 

In Figure 12, we show radial profiles of the gas metallic- 
ity of the SI galaxy cluster, and in Figure 13 we illustrate the 
corresponding mass- weighted gas metallicity maps. When 
additional feedback mechanisms are taken into account, the 
amount of metals in the hot gas component is increased with 
respect to runs with cooling and star formation only (con- 
tinuous blue line on Figure 12 and left panel of Figure 13). 
Also, the metallicity distribution becomes less lumpy be- 
cause the metals are more efficiently transported out of the 
dense regions, increasing the fraction of enriched gas. It can 
be noticed that already bubble heating without winds (red 
dot-dashed line in Figure 12, and middle panel of Figure 13) 
is capable of producing a more homogeneous metallicity dis- 
tribution, while the AGN heating coupled with the galactic 
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winds (green dashed line in Figure 12 and right panel of 
Figure 13) slightly decreases the radial metallicity gradient 
and makes the spreading of the metals even more efficient. 

Nevertheless, the total amount of metals in the three 
different ICM phases - hot ICM gas, cold star-forming gas, 
and stars - remains very similar in all runs, implying that 
there are no substantial metal-enriched gas outflows from 
the galaxy cluster itself. The situation can be rather differ- 
ent if winds with similar velocities are present in less massive 
systems, as demonstrated by Springel & Hernquist (2003). 
When the galactic wind velocities become comparable to 
the escape velocity from the system in consideration, then 
the winds may lead to gaseous outflows, eventually pollut- 
ing the surrounding medium with metals produced by clus- 
ter stars. The heating provided by bubbles might help the 
metal spreading even more. Considering our "BHAR model" 
at early times where it is very efficient and when the cor- 
responding halo mass is considerably smaller, there is evi- 
dence that not only many metals are transfered into the "hot 
phase" , but also that metal enriched gas is pushed towards 
the cluster outskirts, and in part beyond the virial radius. 
However, this effect is transitory, because Ehub decreases 
with time, and more importantly, because the forming clus- 
ter is so massive that it behaves effectively like a closed box 
at late times. Consequently, at z = 0, the distribution of 
metals looks quite similar to the "Magorrian model" dis- 
cussed above. Finally, note that even though bubble heat- 
ing improves metal spreading in the ICM, at the same time 
it does not disrupt the central metallicity gradient of the 
galaxy clusters, as can be clearly seen from Figure 12. There- 
fore, our AGN heating prescription is qualitatively in good 
agreement with the metallicity gradients observed in cool 
core clusters (e.g. Bohringer et al., 2002; De Grandi et al., 
2004; Bohringer et al, 2004). 

4.5 Sound waves or merger induced weak shocks? 

While the unsharp mask technique is very useful in detecting 
X-ray cavities and associated sound waves, it is potentially 
easy to confuse these ripples with shock waves stemming 
from merger events. To demonstrate this danger in inter- 
preting observations, we here analyse a specific case where 
the presence of smaller systems passing trough the cluster 
might induce such a misleading conclusion. 

To this end, we computed projected maps of the S2 
galaxy cluster without AGN heating at z = 0.13. At this 
time, two substructures at radial distances of about half the 
virial radius are moving towards the centre. These substruc- 
tures are still visible in the gas density map (upper right 
panel of Figure 14), but almost completely vanish when the 
X-ray luminosity map is computed. Nevertheless, when the 
unsharp masking procedure is applied to the Lx map, a two- 
lobed feature is clearly visible (upper left panel of Figure 14). 
Even though this feature at first glance looks strikingly sim- 
ilar to the ripples produced by bubble heating events, it is 
exclusively a product of the specific spatial geometry of the 
substructures in the cluster, and is also enhanced due to 
projection effects. 

In the mass-weighted temperature map (lower left panel 
of Figure 14), two spherical regions can be noticed that are 
slightly cooler than their surroundings. They correspond to 
the two substructures. Especially around the substructure in 
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Figure 12. Radial profiles of the gas metallicity of the SI galaxy 
cluster. Only the diffuse hot gas component has been used to 
estimate the metallicity, which is given in Solar units. The blue 
solid line is for the run without AGN heating, the red dot-dashed 
line is for bubble injection based on our "Magorrian model" , while 
the green dashed line is also for the "Magorrian model" but with 
additional inclusion of feedback by galactic winds. We can see 
that the feedback processes manage to better spread the metals 
into the diffuse ICM gas, and the increased mixing leads to a less 
clumpy metallicity distribution. 



the lower right part of the map, hotter surrounding gas can 
be seen. To distinguish between a cold front or a shock, we 
can analyse the pressure of the surrounding gas and compute 
the Mach number map§ (lower right panel of Figure 14). A 
clear jump of pressure in the region adjacent to the hot gas 
around the substructure together with the mildly supersonic 
motion of the substructure indicates the presence of weak 
shocks. Note that these shocks cannot be easily identified 
observationally as being caused by infalling substructures. 
Thus, if a careful analysis is not performed, these features 
could be associated mistakingly with sound waves caused 
by AGN bubbles, especially if at the same time some fossil 
but unrelated bubble is detected in the cluster by its radio 
emission. 

4.6 AGN heating in galaxy clusters of different 
mass 

Here we discuss the effect of AGN heating in clusters span- 
ning a wide range in mass, following their cosmological evo- 
lution up to z = 0. Their main properties are summarized 
in Tables 2 and 3. We analyze how the gas properties of 

§ We made a crude estimate of the Mach number by calculating 
for every gas particle its velocity in the galaxy cluster reference 
frame and dividing it by the local sound speed. 
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Figure 13. Mass-weighted gas metallicity maps of the SI galaxy cluster at z = 0. The left panel corresponds to the case with cooling 
and star formation only. The middle panel shows how the metallicity of the hot gas component changes when AGN heating (here the 
"Magorrian model" was assumed) is included, while the right panel illustrates the case where in addition galactic winds with a velocity 
of ~ 480kms — 1 were included, making the metallicity distribution more homogeneous. 



clusters of different mass and at various radii, normalized to 
i?2oo, change due to the AGN-driven bubbles. 

In Figure 15, we show for all clusters in our sample 
at z — the cumulative X-ray luminosity and mean gas 
temperature at four different radii, r = 0.03i?2oo, O.lifeoo, 
0.5i?2oo and i?2oo, marked with different colours and sym- 
bols. In this Lx — T plane, clusters are found at different 
places according to their mass, with less massive systems 
being located in the lower-left part of the panel and more 
massive ones in the upper-right part, as indicated in the fig- 
ure. Also, cluster luminosity and temperature vary system- 
atically with increasing radius r, with values at the virial 
radius occupying the upper- left part of the plot. The arrows 
attached to each cluster model show how Lx and T change 
when bubble heating is present. Clearly, there is a system- 
atic trend of decreasing X-ray luminosity for all clusters 
and at all considered radii. This effect is more pronounced 
in the most inner regions and it is important both for mas- 
sive clusters and smaller systems. Thus, in the cosmological 
simulations the bubble heating efficiency is not as clearly 
related to the mass of the host cluster as is the case for 
the isolated halos. From Figure 15, it can be seen that the 
two low mass clusters exhibit prominent imprints caused by 
AGN activity, highlighting that the bubble heating is more 
complex in cosmological simulations due to the hierarchical 
merging histories of clusters. 

The gas temperature does not follow an equally clear 
trend as Lx for all the clusters in our sample, but in most 
cases, and especially for central clusters regions, it is boosted 
towards higher values, implying that bubble injection leads 
to en effective heating of the ICM. These trends in Lx and 
T confirm our previous findings for the S1/S2 clusters (pre- 
sented in Figure 9) , and show that they are general features 
of our AGN heating model. Interestingly, recent observa- 
tional work by Croston et al. (2005) indicates that radio-loud 
elliptical-dominated groups have Lx — T scaling relations 
systematically different from those of elliptical-dominated 
radio- quiet groups. Their observed trends in the Lx — T 
plane are qualitatively similar with what we find, show- 
ing that for a given X-ray luminosity radio- loud groups 
have systematically higher gas temperature values. Never- 



theless, since they are probing smaller mass systems, sim- 
ulated galaxy groups with a matching range in mass are 
needed in order to make a more detailed comparison. Re- 
cent observational works (McNamara et al., 2005; Nulsen 
et al., 2005a) revealed the presence of very powerful out- 
bursts due to the AGN activity, which are also accompanied 
with large-scale shocks. These clusters are located above the 
mean Lx — T relation, probably because we are witness- 
ing the very early stages of bubble feedback, where pressure 
equilibrium with the surrounding ICM has not yet been es- 
tablished. 



5 DISCUSSION AND CONCLUSIONS 

In this work, we discussed a simulation model for AGN heat- 
ing in the form of hot, buoyant bubbles, which are inflated 
by active phases of supermassive black holes at the centres of 
massive halos. The motivation for such a mode of feedback 
stems from the rich phenomenology of X-ray cavities and 
radio bubbles observed in clusters of galaxies, and the sug- 
gestion that this AGN activity may represent the solution 
of the 'cooling flow puzzle' posed by clusters of galaxies. 

Several previous studies in the literature (e.g. Chura- 
zov et al., 2001; Quilis et al., 2001; Briiggen, 2003; Hoeft 

6 Briiggen, 2004; Dalla Vecchia et al., 2004) have analysed 
bubble feedback in isolated galaxy clusters using hydrody- 
namical mesh codes. We here present the first implementa- 
tion of this feedback in a SPH code, so an important goal was 
to test whether our results are consistent with these previous 
studies based on very different hydrodynamical techniques. 
Reassuringly, this is the case, both qualitatively and quanti- 
tatively. In particular, for galaxy clusters of similar masses 
as considered by Quilis et al. (2001); Dalla Vecchia et al. 
(2004), and with bubbles parametrised in an analogous way, 
we find changes induced by AGN-heating in gas properties, 
e.g. density and temperature radial profiles, central mass 
deposition rates, that are in excellent agreement. Also, the 
morphology of the bubbles and their time evolution are very 
similar. This is important because it implies that the SPH 
technique, which is more easily applicable to full blown cos- 
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Figure 14. Projected maps of different gas properties of the S2 galaxy cluster at z = 0.13, where a merger with two smaller subclumps 
is in progress. AGN heating has not been included in this simulation. The upper left panel shows the unsharped masked image of the 
X-ray luminosity on a scale of 160 h~ 1 kpc, and the upper right panel illustrates the gas density map. On the lower left panel we show a 
map of mass-weighted cluster gas temperature, while the lower right panel gives a crude estimate of the mass-weighted Mach number. 



mological simulations of cluster formation, can be reliably 
used to study bubble feedback. 

In our simulations, we considered both, isolated halos of 
different mass, and cosmological simulations of the ACDM 
model that follow galaxy cluster assembly from high red- 
shift to the present. The isolated simulations served as a 
laboratory to study the dynamics of bubbles in detail. By 
considering halos with masses ranging from 10 12 /i _1 M© to 
10 15 /i -1 M©, they also allowed us to gain some insight in 
how the coupling of radiative cooling to the AGN heat- 
ing varies with cluster mass. An important conclusion from 
these experiments is that for systems of mass lower than 
~ 10 13 /i _1 M© bubbles with reasonable energy content are 
not capable of preventing excessive gas cooling. If one never- 
theless allows for very large energy in these systems, a stable 
suppression of the cooling requires a delicate fine tuning of 



the bubbles. However, the inefficiency of bubble heating in 
lower mass systems is also caused in part by our injection 
prescription, which is based on global properties of the host 
galaxy cluster without accounting for the actual amount of 
cooling gas present in the very centre. 

In our cosmological simulations, we find that bubble 
injection can substantially affect galaxy cluster properties, 
especially in massive, relaxed clusters and at late cosmolog- 
ical times. Central cluster gas is efficiently heated, and thus 
both the amount of cold baryons and the star formation in 
the central cD galaxy is reduced. Also, an excessive mass 
deposition rate by cooling flows is prevented. AGN-driven 
bubbles not only modify the properties of the most central 
cluster parts, but alter the whole inner region of massive 
clusters, out to radii of order ~ 300/i _1 kpc, where the gas 
density is reduced and the temperature is increased. As a re- 



© 0000 RAS, MNRAS 000, 000-000 



20 D. Sijacki et al. 



10 45 
10 44 


1 1 , i 
- ^ t 

_ / 4 


i 


} * 


1 1 

V 


1 u 


\ $ 




\ 




10 42 






/i x 


\ 


10 41 


\ 






+ 0.03R 200 _ 
* 0.10R 200 
o 0.50R 200 


10 40 


iii 


i 




a 1.00R 200 

, 1 



10 7 10 ! 
T(r) [K] 



Figure 15. Cumulative X-ray luminosity of seven galaxy clusters 
as a function of their mass-weighted gas temperature. We give re- 
sults for redshift z = 0, and at four different radii. Different sym- 
bols and colours denote estimates of Lx and T at different radii 
normalized to R.200'- blue crosses correspond to 0.03i?200, green 
stars to 0.1i?200, yellow diamonds to 0.5i?200, an d red trian- 
gles to i?200- The arrows indicate how the cluster luminosity and 
temperature change when AGN heating is included, i.e. system- 
atically decreasing Lx for all the considered radii. More massive 
clusters are located in the upper-right part of the panel, while less 
massive systems are found in the lower-left corner of the figure, 
as indicated. 



suit, the X-ray luminosity is considerably reduced, while the 
gas entropy exhibits a flat core in the central region. These 
trends are all in the direction required to reconcile hydro- 
dynamical simulations of cluster formation in the ACDM 
model with observations of real galaxy clusters. 

We tried several variants of our bubble model in order 
to explore the dependence of obtained results on the de- 
tailed assumptions made about how the energy is released 
in the bubbles. In particular, we compared an instantaneous 
injection of energy into the bubbles with a scheme where 
the energy is released over a certain period of time (from 
5 x 10 7 yrs to 5 x 10 8 yrs), we imposed that the bubble parti- 
cles should not cool during a given time interval (e.g. ~ 10 8 
yrs) , we tried different spatial patterns for the bubble place- 
ment, and we also varied the initial epoch where our AGN 
heating started (from z = 6 to z = 3) . None of these changes 
was really capable of modifying our results considerably. 
However, it appears that our findings are much more sen- 
sitive to the adopted model for the time evolution of the 
bubble energy. We explicitly demonstrated this by changing 
the rate at which the energy is released with time, under 
the constraint that the total energy injected from z = 3 to 
z = was kept constant. When we linked the bubble en- 
ergy content in this way to a model for the BH accretion 
rate, the fraction of cold gas and the star formation rate can 



be noticeably reduced even at early times, but this is com- 
pensated by a reduced efficiency of bubble heating at late 
times, such that cooling flows are not suppressed sufficiently 
at z = 0. The "Magorrian model", where the feedback oc- 
curs primarily at low redshift fares better in this respect, 
and gives therefore a better match to the properties of ob- 
served rich clusters of galaxies. In addition, we explored yet 
another scaling between the energy of the bubbles and the 
mass of the host galaxy cluster, namely Ehub oc M|qq , which 
can be motivated by the observed Mbh — cr relation as well. 
Also, Ferrarese & Ford (2005) pointed out a relationship 
between the mass of the black hole and that of the host- 
ing dark matter halo, in the form of Mbh oc M^ 5 . Thus, 
if the energy content of the bubbles is a linear function of 
the black hole mass, the above scaling is obtained. One also 
arrives at Ehub oc M^qq if one assumes that the energy in 
the bubbles is some small fraction of the total thermal clus- 
ter energy, which roughly scales as M|qq . In performing the 
analysis with the modified slope, we fixed the normalization 

5/3 

of the Ehub oc M 2 q relation at redshift z = to be equal to 
the value in our ordinary "Magorrian model". The results 
of this analysis showed that the change in the slope from 
4/3 to 5/3 produces qualitatively very similar results, and 
hence it follows that the properties of our simulated galaxy 
clusters are not very sensitive to such a modest change of 
the AGN heating prescription. 

In the newly emerging picture for the joint evolution 
of galaxies and supermassive black holes, the interplay be- 
tween AGN and their host galaxies may be composed of 
two modes. One mode is caused by the quiescent accretion 
of intracluster gas onto the central BH, provoking periodic 
AGN activity which manifests itself in jets and radio bub- 
bles. This mode plays a more important role at late cosmo- 
logical epochs, in massive and more relaxed systems, and it 
is often referred to as a "radio- mode" (e.g. Croton et al., 
2005). The other mode occurs in merging pairs of galaxies, 
where strong tidal forces efficiently funnel large amounts of 
cold gas towards the nuclei of the merging galaxies, where 
it becomes available for fueling the embedded supermassive 
BHs (Di Matteo et al., 2005; Springel et al., 2005b). The 
associated intense accretion triggers quasar activity, which 
is more frequent at higher redshift due to the larger merger 
rates there. If a small fraction of the bolometric luminosity 
of the quasar couples thermally to the surrounding gas, a 
prominent gas outflow can eventually be created during the 
formation of ellipticals, which then shuts off further accre- 
tion and star formation (Springel et al., 2005a) and estab- 
lishes the Mbh — cr relationship (Di Matteo et al., 2005). 

The latter process, the 'quasar mode', has already been 
explored in direct simulation models of galaxy mergers, but 
not yet in cosmological simulations. It would therefore be ex- 
tremely interesting to couple these two modes of AGN feed- 
back in a unified simulation model for supermassive black 
hole growth, and to carry out cosmological simulations with 
it. In such a model, the energetics of the bubbles and the pe- 
riods of AGN activity can then be made directly dependent 
on the current BH mass and on the local physics of accret- 
ing gas, removing much of the freedom in our present bubble 
models. We will present such a model in forthcoming work. 
In addition, our work suggests that for a more complete 
picture of the ICM dynamics, additional physical processes 
should be incorporated as well. This includes thermal con- 
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duct ion, even though it is probably relevant only in the most 
massive systems. We also suggest that the physical viscos- 
ity expected for the ICM gas should be considered as well, 
since the efficiency of non-local AGN heating by viscous dis- 
sipation of sound waves will depend crucially on this input. 
Finally, radio observations strongly suggest that bubbles are 
prevalently filled with relativist ic particles, which appear as 
radio lobes, and in many cases are coincident with X-ray 
cavities. Therefore, it would be important to address, us- 
ing fully cosmological simulation of cluster formation, the 
role of non-thermal radio plasma in bubbles for heating of 
the ICM. It appears that for some time to come clusters of 
galaxies will remain one of the most interesting places to 
study complex hydrodynamic phenomena in the Universe. 
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